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Abstract 

Hamiltonian perturbation theory is used to analyse the stability of f(R) models. The Hamil- 
tonian equations for the metric and its momentum conjugate are written for f(R) Lagrangian in 
the presence of perfect fluid matter. The perturbations examined are perpendicular to R. As 
perturbations are added to the metric and momentum conjugate to the induced metric instabilities 
are found, depending on the form of f(R). Thus the examination of these instabilities is a way to 
rule out certain f(R) models. 
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I. INTRODUCTION 



The question of dark energy has been at the heart of cosmology since the discovery 
of accelerating expansion of the universe The traditional picture of general relativity 
with ordinary relativistic or non-relativistic matter in homogeneous and isotropic universe 
meets severe problems when accommodating it to current cosmological observations. The 
conflicting observational evidence comes mainly^ from supernova light curves [l, 0], CMB 



anisotropies 



and large scale structures 



a. 



This has lead to several suggested 



remedies. Perhaps the most popular way is to add some non-conventional matter to the 
universe. Among these the simplest possibility is no doubt to use the cosmological constant. 
A review of the subject can be found in Q. In any case, the key aspect is the negative 
pressure of the new matter which boosts the expansion of the universe. Other considerations 
include more general distribution of matter, i.e. non-homogeneous or non-isotropic universe 
(see e.g. 

Besides these two, a lot of effort has been put into studies on generalizations and modifi- 



cations of General Relativity. For example metric-affine theories (see e.g. [9(), scalar-tensor 
theory (see e.g. fiol. flit), brane-world gravity (see e.g. 12]) and more general Lagrangians 
have been considered. In the present paper we are especially interested in f(R) gravity 
mode 
R 
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s in which the Einstein-Hilbert action is replaced by a function of the curvature scalar 
I5I If], 17]. None of these modifications is free of problems and this is indeed 
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the case of f(R) gravity as wel 



2fl 
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. As for any model, the cosmological observations issue 
191 ]) as do the observations in the solar system (see e.g. 
24|). The opinions are still divided on the viability of f(R) theories of grav- 



some constraints (see e.g. [18|, 



25, 



261 ] ) as well as sceptical ones (see 



ity There are numerous approving studies (see e.g. 
e.g. 0,0)- 

As the actual universe is not homogeneous and isotropic but contains local perturbations, 



additional challenges for f(R) theories emerge from stability analysis [29|, |30|, 1311 ] . An ac- 
ceptable cosmological model has to be stable against perturbations in the metric and the 
mass distribution. However, stability analysis is customarily done only in the direction of 
R, i.e. only curvature perturbations are considered. This is motivated in particular in the 
case of General Relativity, where the relation between space-time curvature and the matter 
density is a simple one: the trace of Einstein equations imply R oc p. This in turn implicates 
a simple and direct relation between the perturbations in matter and curvature. This is not 
the only possibility. In a f{R) model the model the relation is more complicated due to 
appearance of function f(R) and higher derivative terms in the field equations. The phase 
space is considerably larger and metrics corresponding a given matter distribution ambigu- 
ous. The physical acceptability, however, of a model requires general stability; also stability 
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against perturbations which keep curvature constant, perpendicular to R. 

The Hamiltonian formulation of general relativity has been around since the work of 



Arnowitt, Deser and Misner 
of Ashtekar 
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Hamiltonian formulation has also surfaced in the works 



33| . The first papers on the subject often neglected the boundary terms, 



however, later works have clarified these details (e.g. 34], [35|). Hamiltonian formulation has 
not received too much interest in contemporary papers. In particular and to our knowledge 
the use of Hamiltonian formulation on perturbations of f(R) theories has not been studied 
so far. The main interest has been in specific choices for the function f(R). 

In the present paper we look into perturbations using Hamiltonian formalism of f(R) 
theories. While the technique has not yet been applied to general f(R) theories with per- 
turbations it is a useful tool in studying the stability of f(R) models: with it is simple to 
study perturbations perpendicular to R. As in classical mechanics the Hamiltonian is writ- 
ten as a functional of the fields and their canonical momenta. However, in a geometric theory 
like General Relativity and f(R) theories, some complications appear due to constraints be- 
tween field components. The two main aspects of the canonical Hamiltonian formalism are 
that the field equations are of the first order in the time derivatives and that time is distin- 
guished from other coordinates. For writing the Hamiltonian equations we must thus foliate 
the region of space-time with space-like hypersurfaces. Finally the resulting field equations 
for the perturbations are then analyzed for instabilities. The conventions and details of the 
formalism can be found in 



36|. 



The paper is organized as follows. In section [TT] we write the Hamiltonian field equations. 
The 3+1 decomposition and foliation of the space-time are also presented. The first order 
perturbations are added to the system in section IHIl We also take a look at second order 
perturbations in section HVl In section |V] we summarize our results. 



II. HAMILTONIAN FORMULATION 



In this section we mainly follow the treatment of 36j . Another overview of Hamiltonian 
(and Lagrangian) formulation of general relativity can be found in 37|. As f(R) theories of 
gravity can be written as scalar tensor theories 



S — Sq + s 



20| | . we start by writing the f(R) action as 
]r I f{R)V=g d 4 x + S s + S M 

D7T Jy 

^ J v V=g{f(v) + f(<p)R - fi!P)<p))*x + s s + s M , (i) 



1 The ideas were first seen in the long out of print Gravitation: an introduction to current research. The 
authors have later on released the article on ArXiv as cited. 
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where V is a volume of space-time, S s is a surface term which we will cover later on, Sm is 
the matter term and (p is the auxiliary scalar field. Variation with respect to (p would lead 
to equation f"(<p)(R — ip) =0. Therefore R = <p unless f"(<p) = 0. If f"(p) = the two 
forms of action ([I]) are trivially equal. From here on we make the assumption f"(<p) = 
unless stated otherwise. 

For the purposes of writing the action in terms of the Hamiltonian a 3+1 decomposition 
is needed. We foliate the space-time volume V with space-like hypersurfaces S t of constant 
time. We shall use the Greek alphabet for space-time and the Latin alphabet for space. By 
decomposing the line element to (three) scalar, (three) vector and (three) tensor parts as 

ds 2 = -N 2 dt 2 + h ab {dy a + N a dt)(dy b + N b dt), (2) 

one defines the lapse N, the shift N a , and the induced metric h ab = g^ u e^e b on the hy- 
persurface S t . By using these introduced quantities the invariant volume element reads 
y/—g d 4 x = Nyh dtd 3 x. The Ricci scalar can be written in terms of the extrinsic curvature 
K a b, n a the unit normal to the boundary dV and R, the induced Ricci scalar of the (three 
dimensional) metric h ab as 

R = R + K ab K ab — K 2 — 2V a (V /3 n a n /3 - n a V/). (3) 

Extrinsic curvature is the measure of shrinkage and deformation of an object upon being 
moved a unit interval of proper time into the enveloping space-time. It can be written as a 
function of the induced metric, the lapse and the shift, which appear to be the fields we are 
finally interested in, 

K ab = — (h ab -N a]b - N b]a ) . (4) 

The surface terms of the action ([1]) are not of special interest in this paper. However, it is 
not trivial that these parts do not affect the results. Generally, the surface term must be 



added to the action in order to avoid the need for further boundary conditions (e.g. |38l| ) . 

By choosing the space-time volume V so that, its boundary can be written as a union 
of two space-like hypersurfaces S i2 , — E 4l with normals pointing outwards and a time-like 
hypersurface £>, i.e. dV = £f 2 U (— S tl ) U B. The surface term reads 

Ss— I tf(p)K\h\ x l 2 Sy-Sv 

= iri I f^) K ^ d "y ~ I f\v)KVh d 3 y + [ f'(<p)}CV^ d 3 y) - S , (5) 

where e = n a n a . Here Sq = ^- § dv eK \/\h\ d 3 y is a non-dynamical subtraction term the 
purpose of which is to prevent the integral from diverging in the limit when the spatial 
boundary S t is pushed to the infinity. The constant K is the extrinsic curvature of the 
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boundary dV embedded in flat space-time. In the last term 7 is the induced metric on B 
and K, is the extrinsic curvature scalar of E. However, this is not the only term contributing 
to the surface part. The term f'((p)R from ([!]) produces surface and volume terms, namely 



/ f(tp)R^g~d 4 x = l % dt t f'{<p){R + K ab K ab - K 2 )N^ d 3 y 

JV Jti JT, t 

-2 <f f{ip){V p n a n p -n a V/) dZ a . (6) 

JdV 

When combining these surface contributions the first two terms in ([5]) are eliminated. The 
only surface term left from is 

-2 / f{^){Vpn a n p -n a V^)dT, a = -2 [ f(ip)(V fi n a )n^r a V=ld 3 y 

JB JB 

= 2 / f{ V ){Vpr a )n^n a ^ d 3 y, (7) 



JB 

where r a is the perpendicular unit vector of the boundary St of S t , i.e. r a r a = 1 and 
r a n a = 0. Summing the remaining surface terms we obtain 

S s = 2 I (k - ko)fXv)Nyfc d 2 9. (8) 
J s t 

We have also introduced the induced metric on the boundary St, ctab = habt a A e h B and o is 
its trace. The extrinsic curvature of S t embedded in E t is kAB^A 6 ^ 7 b r a, k is its trace and 
similarly k with the embedding in flat space. The constant k comes from the subtraction 
term. 

We now have the surface part of the action ready for construction of the Hamiltonian. We 
shall see later on that the surface term ([8]) is indeed cancelled in the process of calculating 
ihe field equations. Many of the technical details were omitted and we refer the reader to 



36( | . The generalization to f(R) is easy. 

In the Hamiltonian formulation field equations are found for fields and their momentum 
conjugates. Here the fields are h a b, N, N a and <p. It turns out that in the case of f'(R) 
gravity we need only the momentum conjugate to the induced metric h a b. This can be 
written using the extrinsic curvature 

dh ah oK cd 167r 

For evaluating dK c d/dh a b the extrinsic curvature was written as a function of the induced 
metric given in the formula (jlj). 

For writing the Hamiltonian density Ti. = p ab h a b — ^/—gC we still need the volume part 
of the action. We write the gravitation part of the action without the surface part (which 
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we include later on) as 



SgV = ^St ^ L ^ + f '^ KabRab + R-K 2 )- f'(<p)<p] NVh d 3 x] . (10) 

After some manipulations the volume part of the Hamiltonian density can be cast to the 
form 

H G =P ab h ab -V=jj£ G 



16tt 



{/V) [K ab K ah -K 2 -R + ^\- f(<p)} 



Ho =±. {NVKfrm - m - RfM] + 7 ^ ) {p^ b - ^ )- (13) 

- 2 ^(S),> 3 - (i4) 



+ [(JT* - Kh« b )N a ] |6 - - ^ a6 ),6iVa. (11) 

To express the Hamiltonian density in term of adequate variables, i. e. induced metric h ao 
and its conjugate momentum p ab , we need to rewrite the extrinsic curvature. By inverting 
© we get 

VhK ab f'( V ) = 16n(p ab - X -ph ab ) = p ab - X -ph ab . (12) 

Using this equation the Hamiltonian can be written as a function of the momentum conju- 
gate. Now the volume part of the gravitational Hamiltonian is obtained by integrating Tic 
over the hypersurface S 4 

Jib 

Similarly we get the surface part of the gravitational Hamiltonian by taking the appropriate 
terms and integrating over the hypersurface St: 

' • (h \ I h - L;,) _ 

Vh 

The latter term is produced by applying Gauss theorem to the middle term of ffTTT) when 
integrating over the density. We obtain the field equations by varying the action with respect 
to N, N a , h a },, p a b and ip. As we can see in the action the only time derivatives are those of 
the induced metric. Thus the only dynamic field is h ab and the only momentum conjugate 
needed is p ab . We have the normal boundary conditions for the variations vanishing on the 
boundary 

5N = 5N a = 5h ab = 5<p = 0. (16) 

The full Hamiltonian H includes both surface and volume parts as well as a matter part 
Sm- Since we can write variation of the action as 

^ dt\ f (p ab 5h ab + h ab 5p ab )d 3 y - SH] (17) 



87T 



N (k - ko) - f(^d 2 9. (15) 



5S 
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the Hamiltonian equations are of the form 

• &Hg .06 _&Hg SSm ^c_ n &Hg_S_Sm dHG_ ( , n «x 
ab ~ dp ' P ~ dh + dh' dN a ~ ' dN ~ dN' dip ~ [ ' 

To simplify the field equations, we can choose the foliation to be such that N a = and 
hence h ab = g ab , when the effects of the surface terms vanish. This choice removes one field 
equation, that of N a , and the other ones are much simplified. After tedious calculations we 
end up with equations 

IfaVhp =(R + K 2 - K ab K ab - (p)f'((p)Vh + f((p)Vh, (19b) 
Vhf'(<p) 2 

<p — R=- 5-, (19d) 

where G afc = _R af> — ^Rh ab . For technical details we refer the reader to [3(5!] which can 
straightforwardly be generalized to the f(R) case. Note, however, that in the derivation of 
the equation ( 119dl) further use is made of the assumption f"(R) ^ 0. Otherwise, we would 
get a trivial equality. As can be seen in equations (I19al) and (|19bl) we have also added matter 

SS m Jg^ 5g 00 r- , , 

1W = -f t «'Jn = < 20) 

= _^/l Tab Sg^ = Ny/h phab 
5h ab 2 <ta a6 2 1 J 

which is of the perfect fluid form. 

Even though we assumed from the start that f"(R) ^ it is worthwhile to take look at 
the case of Einstein-Hilbert Lagrangian. If in Eq. we choose f(R) = R, the equality 
is trivial, and only a variation of a constant resulting in a trivial field equation. As the 
assumption of f"(R) ^ is needed in the field equations only in (119dj) the equations would 
stay the same except for this one equation which is irrelevant. From equation f!19bl) we get 
the familiar Friedmann equation for the background 

« 2 = (22) 

oCt 

in a matter dominated universe (p = a~ 3 po)- Here TC = a f (i])/a(r]) is the conformal Hubble 
parameter. We will need this background result later on when we insert the asymptotic 
background solution into the equations. Namely, we can solve a'(i]) from this equation. 
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III. FIRST ORDER PERTURBATIONS 



In this section we add first order perturbations to the metric and the momentum con- 
jugate. In general relativity the trace equation connects curvature and matter density (for 
fixed equation of state p = p(p)) by a simple relation R = n(p — 3P), where k = 8nG. 
The perturbations would be connected correspondingly: 6R = n(6p — 36 P). As the trace 
equation in f(R) gravity is f'(R)R — 2f(R) + 3Df'(R) = n(p — 3P) there are more freedom 
in metrics that produce a given mass configuration. Indeed, the relation between curvature 
and matter distribution is no more an algebraic one, but defined by an differential equation. 
Thus the phase space of metrics is larger and there are perturbations keeping R and thus 
p fixed. This is manifested by the statement that Birkhoff 's theorem^] is no more valid in 



the traditional form in f(R) theories 421 . |43j. Since there are number of studies of the 
perturbations along R (e.g. [29() we are now interested in the opposite and do not introduce 
perturbations to matter but perturbations perpendicular to R only, i.e. 6R = 6p = 0. 

We may add the most general first order perturbations to the metric. These include 
scalar, vector and tensor perturbations. In light of the recent observations and for simplicity 
we examine the case of spatially flat FRW metric. The perturbations in first order can now 
be written as 



44] 



#00 = 9oo - 2a 2 $, (23a) 
9oa = 9oa + a 2 (d a u + u a ), (23b) 
9ab = 9ab + a 2 {- 2^6 ab + V ab x + d a Xb + d b Xa + Xab) , (23c) 

where tilde denotes the background part and the vectors u a and x a are transverse (i.e. 
d a uj a = 0, d a Xa = 0) and Xab is trace free and symmetric tensor (i.e. d a Xab = 0, Xa — 0)- 
Comparing the elements in (123]) and the line element (|2J) to find the perturbations in the 
first order for lapse, shift and the induced metric we obtain 

jV = A + a$, (24a) 

N a = N a + a 2 (d a u + u a ) = N a + a 2 u a , (24b) 
h a b = h a b + a 2 (- 2^6 ab + V a V b x + d a Xb + d b x a + Xab) 

= h ab + a 2 (- 2^6 ah + Xab) ■ (24c) 

The standard practice of splitting the perturbations into scalar, vector and tensor parts 



451 ] is motivated by the reason that in a linear theory these modes decouple. Moreover 



391 ]. However, two years earlier a less known 



2 Birkhoff proved the so called Birkhoff theorem in 1923 
Norwegian physicist Jebsen presented the idea in [40j . The history of the theorem is examinee! in 



4l|. 



s 



each of them has a clear physical interpretation 46]. The first order vector perturbations 
are not generated in the presence of scalar perturbations and dissipate over time. Tensor 
perturbations cause gravitational waves which do not couple to first order scalar perturba- 
tions. Therefore, we may omit vector and tensor perturbations in the first order case and 
assume uj a — 0, Xab = 0. We can further simplify the metric for our purposes by choosing 
an appropriate gauge. We choose to use the Poisson gauge [44| which is a generalization of 
the much used longitudinal gauge. The gauge conditions are 

V • x = 0, (25a) 

V • Cj = 0. (25b) 

Since uo a and Xa are transverse vectors and Xab is a symmetric, transverse and trace-free 
tensor we have u = x = Xa — 0. Along with the physical meaning of the perturbations 
discussed above the perturbed metric simplifies to 

N = N + a$, (26a) 
N a = N a , (26b) 
Kb = Kb - 2a 2 ^5 ab . (26c) 

As the dynamical components of the metric are coupled to their momentum conjugates, 
we are to add perturbations also to the conjugates. Only the induced metric h ab has a 
conjugate p a b, and hence for perturbed one we write 

Pab = Pab + ®5 ab (27) 

having same structure as fl2"6"b). 

In the following we work mostly, unless otherwise stated, in conformal time instead of 
standard coordinate time. So we have ds 2 = —a(rj) 2 dri 2 + a(r]) 2 5 a bdx a dx b , where the confor- 
mal time t] is related to standard coordinate time by drj = a~ x dt. Prime denotes derivatives 
with respect to the conformal time and dot denotes derivatives with respect to the coordi- 
nate time. This choice of background metric corresponds to R = G a t = and \fh, = a 3 . 
Also, we now have p a b = — 2/'(<^)a 3 a'. Since we wrote the f(R) theory using a scalar in ([1]) 
we have <p ~ R. Perturbing (p would produce perturbations parallel to R which we are not 
interested in. 

The Eqs. ( TT9"j) for the chosen background metric and scalar field are now given in a fairly 
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simple form. This reads 



16vrPa 4 = 2(a') 2 f (V?) - « 4 (/(^) - <pf'(<p)) ~ 4a(a"/'(^ + W)), (28a) 

f(^) (6(a') 2 - «V) , 

16vrpa 3 = ^ + a 3 /(<^), (28b) 

a 

6(a') 2 , N 

V = (28c) 

We get only three non-trivial equations as H19c[) produces only a trivial identity. These 
equations, satisfied for any acceptable matter are used to simplify the perturbation equations 
derived later. In the following we assume a matter filled universe with P = and p = po/a 3 . 

By adding the perturbations introduced in (|24p and (|27|) to the equations of motion flTTJ]) 
we get three equations for the large scale perturbations (i.e. space independent perturba- 
tions) 

* = mzm- (29a) 

V a a a°j (ip) J 

$ = 0, (29c) 

where we in (|29bl) the background equations (I28al) and fl28cj) were applied to simplify the 
equation. We immediately notice, that there remains only one dynamic equation while the 
other two are algebraic. The background equation for the induced metric can be used to 
eliminate the second time derivative of the scale parameter. Eq. ( I29bl) is thus written as 

Q'=( 5 --^-)e. (30) 



a a'f'(ip). 

The behaviour of perturbation is clearly dependent on form of the function f(ip) explicitly 
via its derivatives. Moreover, it is found that the time derivative of ^ is zero and therefore 
by equation (I29al) we can write 

6 = Ca 3 a'f{ V ) (31) 

where C is a constant. So, in a universe with growing a(rj) the perturbations in momentum 
conjugate increase. The perturbations of metric tensor, however, behave differently: the 
temporal part vanishes and the spatial perturbations are constant. So, the system leaves 
the linear perturbative regime and ultimately suffers linear instability. 

Although asymptotic analysis is ultimately irrelevant for a linearized unstable system, we 
take a look to some examples to get a better feeling of the evolution. As known, the simple 
function f(R) = R — p 4 /R results asymptotic Einstein-de-Sitter behaviour. Now a(t) = e At , 
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and coordinate and conformal times are related by rj = — ^x~ + c so that a{rf) = A 1 (c—rj) l . 
In the high curvature limit we get 

eM - c 36A4 + /i4 

where C is a constant. The result can also be written more intuitively in coordinate time as 

4At 4 

V ) V 36A 4 ^ ' 

and thus the perturbations increase as time goes to infinity. Here C is another constant. 
Ultimately the first order perturbation theory breaks down; it is not applicable in this case. 
Similar behaviour can seen explicitly for another often used f(R) = R — p?R 2 . 

Even though we have not included perturbations in matter it is worthwhile to check what 
would happen if we did include these perturbations. For a moment we consider p = p + a, 
where a is a perturbation. It turns out that no density perturbations are present, i.e. 
perturbation equation is o = 0. This is not surprising as the matter perturbations are 
coupled to the temporal perturbation of the metric which is also zero. These vanish unless 
(p (which is essentially R) is perturbed, too. 

As we have found, the only dynamical equation is (I29bl) for the momentum conjugate, 
while the two other equations determine, how metric perturbations follow it; they are con- 
straint equations. If these constraints were to be discarded, we would end up with non- 
diagonal perturbations in the metric. Moreover non-existence of temporal perturbations 
is connected with the orthogonality of perturbations to curvature. As it appears that the 
spatial perturbations in the metric do not grow or vanish in time, there is a flat direction of 
phase space, where any spatial first order perturbation is possible and stable. 



IV. SECOND ORDER PERTURBATIONS 



We have now seen that the first order perturbation predicts that f(R) theories suffer 
instability which invalidates first order expansion; equation (I30p reveals that we cannot 
use first order perturbation theory. The next check would be to consider second order 
perturbations, which might give us further understanding of the perturbations involved. We 
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first write the most general form of the metric and the conjugate momentum as 

N = A> + a($ (1) + $ (2) ), (32a) 

2 

N a = N a + a J2( d ^ (r) + u a ] )^ ( 32b ) 

r=l 

K b = h ab + a 2 [ - (2¥« + ¥^)6 ab + £ (V a Va W + dj b r) + d bXa r) + , (32c) 

r=l 

Pafe = Pah + + Q {2) )5 ab , (32d) 



where the upper index z denotes the order of the perturbation. As we have chosen to work 



in the Poisson gauge 45] we have = x^ — Xa = 0. The vector perturbations (j£ 
and x^ab s ^ m remain, however, and some extra attention has to be paid to them. In general 
the scalar, vector and tensor perturbations do not decouple any more in the second order 
perturbation theory. First order vector perturbations contribute to the second order scalar 
perturbations by terms like uo a u a and vice versa. However, first order perturbations do 
not manifest themselves if not present initially. Since we are now interested in to show 
the instability of the system, it is sufficient that some initial condition reveals unstable 
behaviour. In particular we are free to choose initial condition u a (0) = for the vector 
perturbations. First order tensor perturbations can omit them as well. Note, that if we 
were trying to show the stability of the system, the burden of proof would be much heavier: 
we should show, that any choice of initial conditions leads to stable system. 

As mentioned, vector and tensor perturbations in second order cannot be discarded by 
similar arguments. They are strongly affected by first order scalar perturbations. However, 
the second order scalar perturbations are again independent of the tensor and scalar pertur- 
bations of the second order. Therefore, for our purposes, it is sufficient to study only second 
order scalar perturbations, which can be performed rather simply. We write the relevant 
perturbation equations for second order in the same manner as in the previous section. We 
obtain 

^ (2) ' = 0, (33a) 

9,2, -20^( 9,1, ) 2 + 5 ^^»* <2> - < 33b > 
$ (2) = 0. (33c) 

Thus metric perturbation $^ 2 - ) still vanishes and is again constant related to the per- 
turbation of the momentum conjugate G^ 2 -* by (133b ). The perturbation in the momentum 
conjugate is still depending on the form of the f(R). For f(R) = R — /x 4 /R the result is 
the same as in the first order; the perturbation of the temporal part disappear, the spatial 
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part remains constant and the momentum conjugate is has the only dynamical equation. 
It is clear that the instabilities in the first order propagates to the second order as the 
metric perturbations behave in exactly the same way in both first and second order. Thus 
f(R) models may be inherently unstable up to second order when examining perturbations 
perpendicular to R. Because of the similar form of the first and the second order scalar 
perturbations one might conjecture that it is a more general feature of the theory. 



V. DISCUSSION 

Traditionally the stability analysis is performed in the Lagrangian formalism and the 
analysis parallel to R has been carried out before in several papers (e.g. [29]). Many of 
the interesting f(R) models have been found to be inherently unstable in the past 3l|, |47 |. 
However, stability analysis has not yet been used to the full extent as long as the studies con- 
centrate on curvature perturbations only. By using the Hamiltonian instead of Lagrangian 
formulation we examined the large scale cosmological perturbations perpendicular to R with 
non-relativistic matter. These perturbations are fairly easy to examine with the Hamilto- 
nian formulation. The found instabilities are noticeably different to those of previous works 
(e.g. [iilj]). Because of the constraint R = diagonal perturbations of metric and conjugate 
momentum are related to each other. The temporal part of the metric showed up to be con- 
strained by the conjugate momentum one. Moreover, the spatial part of the metric is forced 
to vanish. If these constraints were not satisfied we would have non-trivial perturbations of 
non-diagonal elements of the metrics like goi. 

The perturbations of the momentum conjugate turn out to be the most interesting ones. 
The equation depends explicitly on the form of the function f(R). Some choices of f(R) 
lead clearly unstable cosmological model, but seemingly not all. We have studied some 
well-known f(R) functions and find them unstable. Albeit the physical interpretation of 
the perturbation momentum conjugate is unfortunately not as clear as that of the metric 
perturbations, equation (Q demonstrates the relation between the momentum conjugate 
and the extrinsic curvature. In the 3+1 decomposition the intrinsic curvature R defines how 
the hypersurface is curved whereas the extrinsic curvature defines how each slice is curved 
relative to the enveloping space-time. 

As the perturbations were not well-behaved in this context further studies would be 
relevant in order to find the limits of these constraints. Fruitful directions would likely to be 
investigating the effects of more other types of matter. Also, it would be prudent to examine 
the case where the metric can include shift (i.e. N a 7^ 0). It is clear from the form of (flT|) 
that such a generalization would affect the following equations of motion deeply as the last 
term would be non-zero. This is understandable as the metric would now include spatio- 
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temporal elements. It is also possible to study more general theories with the Lagrangian 
depending also on for example R fJiU R f11 ' or Gauss-Bonnet term. 

It appears that with Hamiltonian formulation of perturbations can be used to constrain 
the spectrum of cosmologically acceptable f(R) theories. While there are several physi- 
cal arguments to judge the f(R) theories like cosmological observations and solar system 
behaviour, stability analysis is one important tool to rule out ill-behaved models out of nu- 
merous possible modified theories of gravity. With continued studies it is possible to find 
the ones best describing the observed behaviour of the universe. 
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